// ProjetoConcorrentes2013.cpp : Defines the entry point for the console application.
//
#include <stdio.h>
#include <math.h>

#define PDF "%9.4f"
#define check(X) \
{ \
        cudaError_t cerr = X; \
        if (cerr != cudaSuccess){ \
                fprintf(stderr, "GPUassert: %s at line %d.\n", cudaGetErrorString(cerr), __LINE__); \
                abort(); \
        } \
}

// aloca mmemoria para vetor de tamanho n
void aloca_vetor(float **V, unsigned int n)
{
    *V = (float *) malloc(n * sizeof(float));
}

// imprime elementos do vetor
void imprime_vetor(float *V, unsigned int n, const char *id)
{
    unsigned int i;

    printf("\n [%s] =\t", id);

    for (i = 0; i < n; i++)
    {
        printf(PDF, V[i]);
    }
    printf("\n");
}

// desaloca memoria utilizada pelo vetor
void desaloca_vetor(float **V)
{
    free(*V);

    V = NULL;
}

// retorna o maior elemento do vetor
__device__ float maxV(float *V, unsigned int n)
{
    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
    unsigned int i;
    float max;

    if (tid == 0)
    {
        max = V[0];

        for(i = 1; i < n; i++)
        {
            if (V[i] > max)
                max = V[i];
        }
    }
    return max;
}

__device__ void calcula_Az(float *A, float *Az, float *D, int n) {
        unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;

        for (int i = 0; i < n; i++) {
                if (tid < n) {
                        Az[i*n + tid] = A[i*n + tid] / A[(n+1)*i];
                        Az[(n+1)*i] = 0.0f;
                        D[i] = A[(n+1)*i];
                }
        }
        __syncthreads();
}

__device__ void calcula_X(float *Az, float *Xa, float *X, float *Bz, int n)
{
        unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;

        float soma = 0;

        if (tid < n) {
                for (int i = 0; i < n; i++) {
                        soma += Az[tid*n + i] * Xa[i];

                }
                X[tid] = Bz[tid] - soma;
        }
        __syncthreads();

        /*
        //reducao do vetor
        soma = 0;
        int i = pow(2, floor(log(static_cast<float>(n))/log(static_cast<float>(2))));
        if (tid >= i && tid < n) {
                V[0] += V[tid];
        }
        i /= 2;
        while (i != 0) {
                if (tid < i) {
                        V[tid] += V[tid + i];
                }
                __syncthreads();
                i /= 2;
        }
        soma = V[0];
        return soma;
        */
}
__global__ void kernel (float *A, float *B, float *Az, float *Bz, float *X, float *Xa, float *difX, float *absX ,int n, float erro, int it_max, int *it, int valor, float *resultado, float *D) {

        unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;
        float mr;

        calcula_Az(A, Az, D, n);

        //calcula Bz: divide B pela diagonal de A
		if (tid < n)
        {
                Bz[tid] = B[tid] / D[tid];
        }
        __syncthreads();

        //copia vetor Bz em Xa para a primeira iteracao (igual especificacao)
        if (tid < n)
        {
                Xa[tid] = Bz[tid];
        }
        __syncthreads();

        //itera para encontrar os valores de X com o erro esperado ou maximo de iteracoes
        *it = 0;
        do {
                calcula_X(Az, Xa, X, Bz, n);

                //coloca os valores absolutos nos vetores difX e absX
                if (tid < n)
                {
                        difX[tid] = fabs(X[tid] - Xa[tid]);
                        absX[tid] = fabs(X[tid]);
                }
                __syncthreads();

                //calcula mr
                mr = maxV(difX, n)/maxV(absX, n);
                //printf("mr = %f\n", mr);      

                //copia X em Xa
                if (tid < n)
                {
                        Xa[tid] = X[tid];
                }
                __syncthreads();

                (*it)++;
		}while((*it < it_max) && (mr > erro));

        //calcula resultado
        if (tid == 0)
        {
                float r = 0;
                for (int i = 0; i < n; i++)
                {
                        r += A[valor*n + i]*X[i];
                }
                *resultado = r;
        }
}


int main(int argc, char* argv[])
{
//      float *Az, *Bz, *X, *Xa, *difX, *absX;
        float *A, *B, resultado = 0, erro;
        float *dev_A, *dev_B, *dev_Az, *dev_Bz, *dev_X, *dev_Xa, *dev_difX, *dev_absX, *dev_D, *dev_resultado;
        int i, j, n, teste, it_max, tam_bloco, it, *dev_it;
        cudaDeviceProp prop;

        cudaGetDeviceProperties(&prop, 0);

        clock_t t = clock();

        scanf("%u", &n);                // dimensao da matriz/vetor
		scanf("%u", &teste);    // equacao de teste
		scanf("%f", &erro);    // valor do erro maximo
		scanf("%d", &it_max);   // numero maximo de iteraçoes

        aloca_vetor(&A, n*n);
        aloca_vetor(&B, n);
/*
        aloca_vetor(&Az, n*n);
        aloca_vetor(&Bz, n);
        aloca_vetor(&X, n);
        aloca_vetor(&Xa, n);
        aloca_vetor(&difX, n);
 		aloca_vetor(&absX, n);
*/
        // leitura da matriz A
		for (i = 0; i < n; i++)
		{
		    for (j = 0; j < n; j++)
		    {
		       scanf("%f", &A[i*n + j]);
		    }
		}
		// leitura do vetor B
		for (i = 0; i < n; i++)
		{
		    scanf("%f", &B[i]);
		}

        //aloca matriz A na GPU
        check(cudaMalloc((void**)&dev_A, n*n*sizeof(float)));
        // copia matriz A para a GPU
		check(cudaMemcpy(dev_A, A, n*n*sizeof(float), cudaMemcpyHostToDevice));

        //aloca matriz Az na GPU
        check(cudaMalloc((void**)&dev_Az, n*n*sizeof(float)));

        // aloca memoria na GPU para o vetor B
		check(cudaMalloc((void **) &dev_B, n * sizeof(float)));
		// copia vetor B para a GPU
		check(cudaMemcpy(dev_B, B, n * sizeof(float), cudaMemcpyHostToDevice));

        // aloca memoria na GPU  para o vetor Bz
    	check(cudaMalloc((void **) &dev_Bz, n * sizeof(float)));

        // aloca memoria na GPU para o vetor X
    	check(cudaMalloc((void **) &dev_X, n * sizeof(float)));

        // aloca memoria na GPU para o vetor Xa
    	check(cudaMalloc((void **) &dev_Xa, n * sizeof(float)));

        // aloca memoria na GPU para o vetor difX
    	check(cudaMalloc((void **) &dev_difX, n * sizeof(float)));
                                                                  
        // aloca memoria na GPU para o vetor absX
    	check(cudaMalloc((void **) &dev_absX, n * sizeof(float)));

        // aloca memoria na GPU para o resultado
        check(cudaMalloc((void**)&dev_resultado, sizeof(float)));

        // aloca memoria na GPU para o resultado
        check(cudaMalloc((void**)&dev_it, sizeof(float)));

        //aloca memoria na GPU para o vetor diagonal (auxiliar)
        check(cudaMalloc((void**)&dev_D, n*sizeof(float)));


        dim3 nthreads(prop.maxThreadsPerBlock);
        tam_bloco = ((n*n % prop.maxThreadsPerBlock != 0) ? (n*n/prop.maxThreadsPerBlock+1) : (n*n/prop.maxThreadsPerBlock));
        dim3 nblocos(tam_bloco);
        /*
        if (n < prop.maxThreadsPerBlock) {
                dim3 nthreads(250);
                tam_bloco = ((n*n % 250 != 0) ? (n*n/250+1) : (n*n/250));
                dim3 nblocos(tam_bloco);
        }
        else {
                dim3 nthreads(prop.maxThreadsPerBlock);
                tam_bloco = ((n*n % prop.maxThreadsPerBlock != 0) ? (n*n/prop.maxThreadsPerBlock+1) : (n*n/prop.maxThreadsPerBlock));
                dim3 nblocos(tam_bloco);
        }
        */
        printf("max threads: %d, nblocos: %d\n", prop.maxThreadsPerBlock, n*n/prop.maxThreadsPerBlock+1);
        kernel<<<nblocos, nthreads>>>(dev_A, dev_B, dev_Az, dev_Bz, dev_X, dev_Xa, dev_difX, dev_absX, n, erro, it_max, dev_it, teste, dev_resultado, dev_D);

        /*
        //DEBUG
        check(cudaMemcpy(Az, dev_Az, n*n*sizeof(float), cudaMemcpyDeviceToHost));
		check(cudaMemcpy(A, dev_A, n*n*sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(B, dev_B, n*sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(Bz, dev_Bz, n*sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(X, dev_X, n*sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(Xa, dev_Xa, n*sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(difX, dev_difX, n*sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(absX, dev_absX, n*sizeof(float), cudaMemcpyDeviceToHost));


        imprime_vetor(A, n*n, "A");
        imprime_vetor(Az, n*n, "Az");
        imprime_vetor(B, n, "B");
        imprime_vetor(X, n, "X");
        imprime_vetor(Bz, n, "Bz");
        imprime_vetor(V, n, "V");
        
        imprime_vetor(Xa, n, "Xa");
        imprime_vetor(difX, n, "difX");
        imprime_vetor(absX, n, "absX");
        
        */


        t = clock() - t;

        check(cudaMemcpy(&resultado, dev_resultado, sizeof(float), cudaMemcpyDeviceToHost));
        check(cudaMemcpy(&it, dev_it, sizeof(float), cudaMemcpyDeviceToHost));

        printf("\n---------------------------------------------------------\n");
        printf("Iterations: %d\n", it);
        printf("RowTest: %u => [%f] =? %f\n", teste, resultado, B[teste]);
        printf("---------------------------------------------------------\n");
        printf("\n Tempo de execucao : %f s \n\n", t / (double) CLOCKS_PER_SEC);

        desaloca_vetor(&A);
        desaloca_vetor(&B);
/*
        desaloca_vetor(&Az);
        desaloca_vetor(&Bz);
        desaloca_vetor(&X);
        desaloca_vetor(&Xa);
        desaloca_vetor(&difX);
        desaloca_vetor(&absX);
*/
        cudaFree(dev_A);
        cudaFree(dev_B);
        cudaFree(dev_Az);
        cudaFree(dev_Bz);
        cudaFree(dev_X);
        cudaFree(dev_Xa);
        cudaFree(dev_difX);
        cudaFree(dev_absX);
        cudaFree(dev_it);
        cudaFree(dev_resultado);
        cudaFree(dev_D);
        //printf("\n\nFIM!\n\n");

        return 0;
}
                                                                                                                             

